Introducción

En esta práctica vamos a tratar con datos de defunciones por enfermedad isquémica en Aragón, separados por municipio. Nuestro objetivo va a ser modelizar el riesgo de muerte a causa de esta enfermedad, y realizar una predicción del riesgo “real” de muerte para cada municipio.

1. Lectura y preparación de los datos

# Librerías
# install.packages("rmdformats") # plantillas para reporte HTML
library(dplyr) # operaciones con tablas
library(sf) # lectura y tratamiento de datos geoespaciales
library(spdep) # lectura y tratamiento de datos geoespaciales
library(R2WinBUGS) # modelos bayesianos, llama a WinBUGS
library(INLA) # modelos bayesianos utilizando INLA
library(ggplot2) # representación gráfica
theme_set(theme_void()) # cambiar estética de todos los plots de ggplot
library(gridExtra) # extensión de ggplot para plotear varias gráficas a la vez
library(leaflet) # representación de mapas interactivos

Comenzamos cargando las librerías necesarias para nuestra práctica y leyendo los archivos de datos:

# Cargamos los datos con info por municipio
load(file.path("datos", "Aragon.Rdata"))
datos <- Aragon.df

# Leemos los datos georeferenciamos y los unimos con los datos cargados
mapa <- read_sf(file.path("datos", "aragon.shp")) %>%
  left_join(datos, by = "CODMUNI") %>%
  mutate(RME = O / E)

Para cada municipio tenemos el número de muertes observadas por enfermedad isquémica en hombres, así como el número esperado de muertes, calculado como el número de muertes que le correspondería a cada municipio dada su población, multiplicando la población del municipio por el número promedio de defunciones por habitante de Aragón. Como una primera aproximación, vamos a representar los observados, los esperados y el riesgo de mortalidad estandarizado en un mapa:

Se puede apreciar en los mapas del número de casos observados y esperados que el número de fallecimientos es muy dependiente de la población de cada municipio, y que si no lo ajustamos de alguna forma lo único que vamos a observar en el mapa es un reflejo de la población de cada municipio. Al dividir observados entre esperados obtenemos la Razón de Mortalidad Estandarizada (RME), que ya es un indicador más preciso del riesgo real de fallecimiento. Un RME mayor a 1 indica un exceso de riesgo, por lo que el mapa en el que representamos el RME es útil para conocer las zonas de mayor riesgo. El problema es que muchos municipios de Aragón poseen muy pocos habitantes, por lo que una desviación muy pequeña entre observados y esperados causa variaciones muy grandes en el RME, teniendo zonas con un riesgo aparente de 0, y zonas con un exceso de riesgo enorme, cuando realmente sólo hay uno o dos fallecimientos de diferencia entre observados y esperados. Es por esto que vamos a ajustar un modelo que nos suavice este riesgo, teniendo el cuenta los casos observados en los municipios vecinos, para tratar de encontrar una estimación del riesgo que se asemeje más a la realidad.

2. Ajuste de modelo de Besag, York y Mollié

Vamos a implementar el modelo de Besag, York y Mollié (BYS), tanto en WinBUGS como en INLA. Este es un modelo jerárquico que combina un modelo poisson gamma con un efecto autorregresivo que confiere un efecto espacial, al tener en cuenta las observaciones en las localidades colindantes.

# Generamos estructura de vecindad
vecinos <- poly2nb(mapa) %>%
  nb2WB()

# Modelo BYS en WinBUGS
set.seed(1234)

modeloBugs <- function(){
  # verosimilitud
  for(i in 1:n){
    O[i] ~ dpois(mu[i])
    log(mu[i]) <- log(E[i]) + m + het[i] + sp[i]
    het[i] ~ dnorm(0, prechet)
    R[i] <- exp(m + het[i] + sp[i]) # RMS
  }
  for (i in 1:n) {
    pR[i] <- step(R[i]-1) # probabilidad de que riesgo superior a 1
  }
  
  # distribuciones a priori
  sp[1:n] ~ car.normal(adj[], w[], num[], precsp)
  m ~ dflat()
  prechet <- pow(sdhet, -2)
  precsp <- pow(sdsp, -2)
  sdhet ~ dunif(0, 10)
  sdsp ~ dunif(0, 10)
  
  # distribución predictiva
  for (i in 1:n){
    yPred[i] ~ dpois(mu[i])
    res[i] <- yPred[i] - O[i]
    pRes[i] <- step(res[i])
    pY[i] <- step(yPred[i])
  }
}

datos<-list(O = mapa$O,
            E = mapa$E,
            n = nrow(mapa),
            adj = vecinos$adj,
            w = vecinos$weights,
            num = vecinos$num)

iniciales<-function(){
  list(m = rnorm(1),
       het= rnorm(nrow(mapa),0,0.1) ,
       sp = rnorm(nrow(mapa),0,0.1),
       sdhet=runif(1,0,1),
       sdsp=runif(1,0,1),
       yPred = rpois(nrow(mapa),mean(mapa$O)))
}
param <- c("R","m", "het","sp","sdhet","sdsp","pR","pRes","yPred","pY")

resultados <- bugs(
  data = datos,
  inits = iniciales,
  parameters = param,
  model = modeloBugs,
  n.burnin = 3000,
  n.thin = 25,
  n.iter = 10000
  # debug=TRUE
)

# Generamos estructura de vecindad y la guardamos en un archivo
# poly2nb(mapa) %>%
#   nb2INLA(file = file.path("datos", "aragon.graph"))

# Modelo BYS en INLA
H <- inla.read.graph(filename = file.path("datos", "aragon.graph"))

S <- U <- seq(1, 729)

mapa <- mapa %>%
  mutate(S = S,
         U = U)

formula <- 
  O ~ 1 + 
  f(S, 
    model = "besag",
    graph = H,
    scale.model = TRUE,
    hyper = list(
      prec = list(prior = "loggamma", param = c(1, 0.001)))
  ) +
  f(U,
    model = "iid",
    hyper = list(prec = list(prior = "loggamma", param = c(1, 0.001)))
  )

modeloINLA <- inla(
  formula = formula,
  family = "poisson",
  data = mapa,
  E = E,
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
  control.predictor = list(compute = TRUE, cdf = c(log(1)))
)

3. Mapa interactivo con librería leaflet

Vamos a representar en un mapa interactivo todos los riesgos y probabilidades estimados, para tener toda la información obtenida a lo largo de la práctica en un único lugar y poder visualizarla y compararla más eficientemente.

Reproducibilidad

La elaboración de este informe ha sido gestionada mediante el uso de git. Puedes descargar todo el código y los datos en este link, o, si eres usuario de git, utilizando el comando git clone https://github.com/Julian-Guillo/Disease-Mapping-Aragon.

Con respecto al código, este ha sido extensamente comentado, los nombres de las variables son representativos de su contenido y siguen la convención camelCase. Para cualquier corrección o duda por favor abrir un pull request o contactar conmigo a julianleioa@gmail.com.